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SUMMARY 

This  report  discusses  the  implementation  of  system  and 
process  noise,  including  sensor  errors  and  atmospheric 
turbulence,  in  the  simulation  of  continuous  models  by 
digital  computers.  The  concept  of  discrete  white  noise  is 
first  introduced  in  a  manner  in  which  it  reduces  to 
continuous  white  noise  as  the  integration  time  interval 
reduces  to  zero.  The  derivation  of  first  and  second  order 
Markov  processes  from  white  noise  is  then  discussed.  The 
discussion  includes  the  consideration  of  a  suitable 
difference  algorithm  to  approximate  differentiation.  The 
treatment  throughout  has  been  aimed  at  providing  the 
reader  with  the  tools  to  implement  a  required  power 
spectral  density  (or  autocorrelation  function)  as  noise  in 
any  discrete  digital  simulation  model. 
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1.  INTRODUCTION 

In  this  Report  we  discuss  the  implementation  of  noise  in  the  discrete 
simulation,  on  a  computer,  of  a  continuous  process.  We  will  develop  the 
concepts  necessary  for  our  purposes,  but  only  in  sufficient  depth  to  preserve 
continuity  of  the  text.  References  are  given  to  books  in  which  a  rigorous 
mathematical  treatment  can  be  found. 

In  the  next  section  we  briefly  examine  continuous  noise,  starting  with  white 
noise  and  then  showing  how  noise  of  any  desired  spectral  density  can  be 
obtained  by  using  white  noise  as  the  input  variable  to  a  suitable  (stochastic) 
differential  equation. 

Section  3  examines  discrete  noise,  using  the  z  transform  formalism,  in  a 
manner  suitable  for  implementation  in  digital  computer  simulation  of  a  syster. 
Particular  care  is  taken  to  develop  a  formalism  which  reduces  to  the 
continuous  case  as  the  discrete  interval  tends  to  zero.  The  z  transform,  as 
normally  defined,  does  not  have  this  property. 

The  discrete  equivalent  of  white  noise  is  discussed,  and  then  the  spectrum  of 
Markov  processes  generated  by  using  white  noise  as  the  forcing  function  to  a 
difference  equation  is  examined.  The  relative  merits  of  various  algorithms 
for  the  discrete  approximation  of  the  continuous  differential  are  considered 
in  this  treatment. 

Approximating  a  continuous  noise  process  by  an  equivalent  discrete  sequence 
introduces  errors,  which  increase  with  the  integration  step  size.  Expressions 
are  derived  relating  the  maximum  step  size  to  the  fractional  error  in  the 
power  spectral  density  and  the  parameters  of  the  system. 

Appendix  I  gives  a  worked  example  evaluating  the  parameters  of  a  given  noise 
spectrum,  and  implementation  in  a  discrete  simulation.  Appendix  II  shows  how 
atmospheric  turbulence  can  be  simulated  by  the  methods  described  in  the 
Roport.  Appendix  III  discusses  the  discrete  simulation  of  the  Wiener  process. 


2.  NOISE  IN  CONTINUOUS  SYSTEMS 
2.1  Methods  for  describing  stochastic  processes 

I f  v  is  some  variable,  eg  acceleration  roll  rate  or  wind  speed,  and  m  is 
the  measurement  of  that  variable,  then  those  are  related  by 


m  =  v  +  n 


(H 


where  n  is  the  measurement  noise  or  turbulence,  represented  by  a  stochastic 
process.  n  is  described  by  the  probability  density  function  p(a,t),  or 
equivalently,  by  the  statistical  moments.  Gaussian  random  processes  are 
completely  described  by  the  first  and  second  moments,  and,  since  most 
random  processes  arc  found  to  be  nearly  Gaussian,  it  is  common  not  to 
exar-ine  any  moment  higher  than  the  second.  The  first  moment,  or  mean,  is 
given  by: 


H  ■  Eln(t)] 


f 


a 


p(a,t)da 


Lt  \ 
T*°°  T 


n(t)dt 


'  -T/  2 


.ao 


(2) 
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where  p  is  the  probability  density  function  for  the  process  with  amplitude 
a  at  time  t.  The  last  equality  follows  from  the  assumption  that  the 
process  is  ergodic.  In  this  memorandum  we  assume  the  noise  is  stationary, 
so  that  p  does  not  depend  on  time. 

The  second  order  moment  can  be  given  in  terms  of  either  the  autocorrelation 
function  or  the  covariance  function. 

/*  00 

C(r)  =  E[n(t  +  t)  n(t)]  =  /  u0p(u,  t;  0,t  +  r)da  $ 

J  _oo 


Lt 

T  oo 


f  T/2 

|  n(t  +  r)  n(t)dt 
J  -T/2 


(3) 


R(r)  =  E[{n(t  +  t)  -Mifn(t)  -  (4) 

where  p(a,t;3,t+t)  is  the  joint  probability  distribution  function.  The 
autocorrelation  and  covariance  functions  are  related  by 


R(t)  =  C(t)  -  v2 


(5) 


In  this  paper  we  will  be  solely  concerned  with  zero  mean  processes  for 
which  V  =  0.  When  dealing  with  a  non-zero-mean  process,  it  is  assumed  that 
the  mean  is  first  removed  from  the  process. 

The  second  moment  can  ^e  equivalently  described  by  the  power  spectral 
density(rof . 1) : 


C(T)e 


-  jcj  r 


dr 


(6) 


which  is  the  Fourier  (or  bilateral  Laplace  transform  with  s  =  ju)  transform 
of  the  autocorrelation  function.  Note  that  we  are  using  a  two-sided 
spectrum.  This  will  he  the  case  throughout.  The  inverse  of  equation  (6) 
is 


C(T) 


1_ 

25T 


.* 

I 


jtn  r 

'h(cJ)  e  dw 


(7) 


Wo  will  now  examino  Pnrsoval  s  theorem,  in  preparation  for  comparing  its 
nxpt oss ions  in  tho  discrete  and  continuous  cases.  Applying  Parsoval's 
theorem  to  Rp  the  Fourier  transform  of  itj.,  where  =  n  for  -T/2  <  t  <  T/2 

and  ts  o  for  ltl>T/2,  gives  (ref  .1) : 


Lt 

T  ♦  oo 


_1 

T 


ni(t)dt 


Lt 

X  +  00 


1 

2sT 


I 

*  _co 


NT(-tu)dto 


(8) 


Thio  theorem  can  bo  regarded  as  a  statement  of  equal  itv  of  energy  in  tho 
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baud  side  of  equation  (8)  is  the  .autocorrelation  at  zero  lag,  ie  the 
variance.  Also,  it  can  be  shown  that 


fNTMNT('u) 


0) 


so  that  equation  (8)  ran  bo  written  as 


oo 

C(0)  =  f  *(“•>)  dw  (10) 

J  _oo 

Equation  (10)  also  follows  from  equation  (7)  by  setting  T  =  0;  however,  we 
have  derived  it  using  Parseval's  theorem  for  comparison  with  the  discrete 
analysis  of  the  next  section.  Note  that  for  continuous  white  noise  both 
sides  of  the  equation  are  infinite  unless  $(w)  =  o. 

2.2  White  noise  and  the  generation  of  Markov  processes 

A  stochastic  process  commonly  used  in  the  study  of  noise  is  Gaussian  white 
noise  which  has  autocorrelation  and  power  spectral  density  defined  by 


Cu(T)  =  Y  6d(t),  <I>u  («)  =  Y 


(U) 


whore  5^  is  the  Dirac  impulse  function,  and  Y  is  the  strength  of  the  noise. 

We  :  ;  j.  *•<  '"or  *-o  this  process  as  continuous  white  noise,  to  distinguish  it 
frv.  :>.  i  :te  white  noise  introduced  in  the  next  section.  The  words 
’continuous1  and  discrete’  are  used  to  describe  the  system  to  which  the 
noise  is  applied. 

If  wo  have  a  measurement  noise  which  is  other  than  white,  this  random 
process  can  be  generated  by  having  it  as  the  dependent  variable  in  a 
differential  equation  whose  forcing  function',  is  white  noise,  and  augmenting 
this  differential  equation  to  the  system  equations (ref .2).  For  example,  if 
x  is  given  by 


x  =  -ux  ♦  au 


(12) 


where  u  is  white  noise  whose  autocorrelation  and  spectral  density  is  given 
by  equation  (11),  then  x  is  also  a  Gaussian  process  whose  autocorrelation 
and  spectrum  are  given  by 


Vr> 


_a.  c'4,lrl  * 
2p 


♦  uJ 


1  * 


(13) 


There  is  an  arbitrariness  in  the  definition  of  a  and  u  in  equation  (12). 
We  shall  remove  this  be  defining  the  white  noise  to  have  unit  power 
spectral  density,  viz  Y  =  1. 
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The  process  defined  by  equation  (12)  is  a  first  order  Markov 
process(ref .2) .  This  noise  is  also  known  as  exponentially  correlated  noise 
from  the  shape  of  its  autocorrelation  function. 


A  second  order  Markov  process 


is  generated  by  the  pair  of  equations: 


y  =  cooX-2fco0y  + 


cu 


(14) 


x  =  y  +  au 


(IS) 


Note  that  the  same  white  noise  source  is  input  to  both  equations. 

The  autocorrelation  function  and  power  spectral  density  of  x  which  result 
from  this  process  are: 


Cx(T)  =~4F^T'cost,  coslvH  -  r  .  Wolrl 


(16) 


♦iM 


a3  to  3  ♦  bl 


oj  4  ♦  2w  I  (2£2  -  l)w  2  ♦  u  S  •  ^ 


(17) 


where,  once  more,  u  is  taken  to  have  unit  power  spectral  density,  and 

b  =  c  +  2ti£u  . 

o 


3.  DISCRETE  SIMULATION  OK  NOISE 
3.1  Methods  for  describing  discrete  random  processes 

In  a  digital  computer  program  which  simulates  a  continuous  system,  the 
states  are  evaluated  only  at  discrete  time  intervals.  For  the  case  of 
uniform  sampling  intervals,  A,  equation  (!)  becomes 


Kv 

*5$ 

i  v* 


4 

il 

m 


mfki)  =  v(kA )  *  iiikA) 


{  U) 


whereby  the  values  of  the  continuous  functions  are  only  sampled  at  set 
times.  The  formal  expression  of  sampling  may  be  written  in  diffetent  ways. 
One  r.caoon  method  (see,  for  example,  reference  3)  is  to  introduce.  Dirac, 
delta  functions  to  represent  the  sampling  process.  However,  this  is 
unsatisfactory  for  computer  simulations,  which  require  numbers  and  not 
impulses.  In  addition,  by  expressing  sampling  in  terms  of  numbers,  as  is 
done  in  equation  (IS),  the  analysis  reduces  to  the  continuous  case  as  the 
sampling  interval  goes  to  zero.  This  is  further  discussed  after 
equat ion  (21). 

Given  the  sampling  process  expressed  by  equation  (181.  the  question  which 
Isas  to  be  answered  is:  uhat  is  tne  relation  between  the  antocorre  5  at  ion 
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function  and  power  spectrum  of  this  discrete  system  and  of  the  continuous 
system  to  which  it  approximates?  We  commence  by  defining  the 
autocorrelation  function  and  spectral  density  for  discrete  time  random 
processes: 

C (lA)  =  Eln(jk  '/{A)  »(kA)) 


Lt  _  1_ 

K  ♦  00  K  +  1 


n(Jk  ♦  l |A)n(kA) 


k=-K/2 


=  A  )  C(*A)z_* 


where  A  is  the  discrete  time  interval  between  samples.  The  power  spectral 
density  is  identified  as  the  two  sided  z  transform! ref . 4 , 5 )  of  the 
autocorrelation  function.  Note  that  the  spectral  density  of  this  sampled 
process  is  only  defined  for  -u/A  <  u  <  tf/A.  The  inverse  of  equation  (20) 
is 


C(*A) 


l- 1 

§  $(2)2  dt 


which  may  be  also  written  as 


by  substituting  7,  -  a 


% 

i 

IT /A 

2X 

1 

J 

JuA 

'  -ir/A 

(21(a)) 


The  factor  A  outside  the  summation  in  equation  (20)  and  outside  the 
integral  in  equation  (21)  has  been  included  to  give  continuity  with  the 
Laplace  transform  as  A*Q,  kA*l.  This  factor  is  not  generally  included  in 
the  definition  of  the  2  transform,  but  is  vital  for  the  treatment  of 
simulation  in  this  test.  We  will  therefore  digress  briefly  on  this  matter 
by  means  of  example.  Consider  the  damped  exponential  defined  by 


x  =2  exp  (-at) 


t  £  0 


vhose  Laplace  transform  is 


t  ■-  0 


X  =  /  jxj  =  (a  +  s)'1 


Taking  the  2  transform  of  x  using  the  A  multiplied  definition 
(equation  (20))  gives 
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*z  =  K|X| 

A  A  *  | x* | 

=  A/ { 1  -  cxp(-(a  +  s)A|),  z  -  csA 


where  x*(t)  =  Ix(t)6p(c-nA) . 

As  A^O.  we  have  X  ,  which  would  not  have  been  the  case  if  the  A  factor 
z  s 

had  not  been  included.  In  general  (eg  in  sampling  systems),  this 
continuity  from  the  z  transform  to  the  bap  lace  transform  is  not  required, 
so  the  factor  A  is  omitted.  However,  we  feel  that  the  z  transform 
including  a  A  factor  is  the  more  nitural  definition.  Confirmation  of  this 
is  suggested  by  considering  the  express  ion  for  the  z  transform  as  a  sum  of 
the  values  of  the  Laplace  transform  over  the  primary  and  complementary 
strips.  Using  the  usual  definition  of  the  z  transform  this  relation  is 


XZ  n  i  )  ys  * 2,rjn/A) 


„CO 


If  the  A  multiplied  expression  for  the  2  transform  is  used,  the  above 
relation  becomes 


y.  v  * 


2*jn/A) 


-OO 


As  AM),  X  (s  +  jn.2*/A)*0  for  all  it  axeopt  n  =*  0  and  wc  ar<*  left  with  the 

s 

simple  relation 


I XJ 


A  *0 


X 


when  we  use  the  A  multiplied  definition  of  the  n  transform.  We  stay 
heuristics! ly  view  the  A  factor  as  moderating  the  impulses  in  the  sampling 
process,  so  that,  as  A'*!!,  we  regain  a  continuous  function.  The  benefits  of 
the  A  multiplied  definition  will  become  manifest  in  Section  3.2  where  si 
enables  us  to  define  a  discrete  equivalent  of  continuous  white  noise.  As 
we  have  said  earlier,  the  A  factor  is  usually  omitted,  since  it  is  of  no 
consequence  when  this  limit  is  not  of  interest. 


Further  insight  to  this  definition  of  the  ?.  transfore  can  he  gained  by 
considering  Pavsevai's  theorem,  which  relates  the  energy  in  the  signal  to 
the  energy  in  the  spectrum.  The  discrete  equivalent  of  Parseval's  theore® 
is 


Lt  A 


K/2 

V  - 

) 

u 

•M2. 


n*  fkA)  = 


Lt 


_  1 
2*K 


.T/& 


-x/L\ 


N*k(c 
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Wo  identify  the  1 1*  f  L  hand  side  with  eorrolction  at  zero  shift  (vide 
equation  (I*)))  and  it  can  be  shown  that 

*U)  =  K^If  NK(z)  NK(z'1)'  z  =  ejWA  (23' 

where  Nv  is  the  z  transform  of  n(kA),  k  =  -K/2  to  K/2. 

Substituting  equatio:.  (23)  into  (22)  yields  equation  (21(a))  with  1=0. 

Contrary  to  the  case  for  continuous  white  noise,  harsoval's  theorem  does 
have  application  in  the  discrete  equivalent  of  white  noise,  as  we  shall 
see . 

3.2  Discrete  white  noise 

Suppose  for  our  noise  in  equation  (18)  we  generate  random  numbers  u(kA) 
with  a  Gaussian  distribution  and  variance  o2.  This  can  be  regarded  as 
discrete  'white'  noise  with  autocorrelation  function  and  power  spec*  ra 1 
density  given  by 


Cu(£A)  =  o26KrU),  $uU)  =  02A  (24) 

where  6  is  the  Kronecker  delta.  We  wish  to  show  that  u(kA)  reduces  to 
Kt* 

continuous  white  noise  defined  by  equation  (11)  as  A'*0(ref.l).  Although 
the  steps  which  follow  are  not  mathematically  rigorous,  they  are  included 
as  an  heuristic  derivation  of  the  required  relations. 

If  we  define  u(t)  as 


uCt)  3  (25> 

whare  the  limit  is  assumed  to  exist  in  some  generalized  sense,  then  we  wish 
to  show  that  u(t)  has  the  properties  of  continuous"  white  noise  as  defined 
hy  equation  (11).  ^ow 


(r )  *  E{  u(t)  u(t  *  f)i 

s  4*G,k4»t,M*T  u^k 

s  ^0,|r-^<  A/2  °5  6Kr(M) 

*  *0  -  W)  -  h(r  ♦  4/2)} 

=  ,l[o!  AUr) 
a*0  D 


where  hit)  denotes  the  unit  step  f-,  sect  son .  Comparing  this  rests;*  with 

equat ion  f 1! ) : 
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n 

II 


(T) 


n,](T) 


nn 


wc  see  thill  u(i)  1ms  the  character  1st. Lcs  of'  continuous  white  noise.,  with 


T  -  A*0  “2A 


(28) 


io  we.  have  shown  that  if  wo  implement  discrete  white  noise  with  the 
autocorrelation  function  given  by  o*6^r(t)  then  this  is  equivalent  to  a 

continuous  white  noise  process  with  autocorrelation  function  Aos6^(x). 


Let  us  continue  this  equivalence  of  continuous  and  discrete  processes  into 
the  frequency  domain.  Continuous  white  noise  has  the  uniform  power 
spectrum  Y  (equation  (11)).  The  discrete  white  noise  has  the  spectral 
density  0*A,  which  extends  from  -tr/A  <  w  <  ir/A.  Using  relation  (28)  it  is 
seen  that  the  amplitude  of  the  power  spectral  density  function  for  the 
discrete  case  becomes  Y,  which  is  identical  to  that  for  the  continuous; 
these  relations  are  shown  in  figure  t.  Note  that  the  discrete  equivalent 
of  Parseval's  theorem  is  obeyed,  so  our  use  of  the  factor,  A,  has  been 
consistent . 


3.3  Discrete  algorithms  for  the  differential  operator 

As  can  be  done  for  continuous  processes,  it  is  possible  to  generete  a 
discrete  random  process  with  a  desired  power  spectral  density  by  a 
difference  equation  with  a  disciete  white  noise  forcing  function.  Since  we 
will  be  wishing  to  relate  the  continuous  and  discrete  processes  thus 
produced,  in  particular,  to  see  how  well  the  power  spectral  density  of  the 
latter  approximates  to  that  of  the  former,  we  will  briefly  digress  and 
examine  the  derivation  of  difference  equations  from  differential  equations. 

The  simplest  discrete  approximation  to  the  derivative  is  the  forward 
difference  defined  by: 


x[(n+l)A]  =  x[nA]  +  x[nA].A 


(29) 


The  difference  equation  thus  obtained  is  unsatisfactory  for  an  integration 
routine,  as  solutions  obtained  using  it  rapidly  diverge.  This  is  because 
the  magnitude  of  the  transfer  function  is  greater  than  the  1/w  value  of 
continuous  integration,  especially  at  the  higher  frequencies.  This  is 
shown  in  figure  2  which  plots  the  ratio  of  the  magnitude  of  the  transfer 
function  for  this  case, 


A  _  juA 

V  2(1  -  cosfl)'  z  °  e 


(30) 


as  a  function  of  angular  frequency. 

All  this  is  well  known,  and  forward  differences  are  rarely  used  in  discrete 
integration  routines.  Commonly  used  integration  algorithms  are  fourth 
order  Runge-Kutta  and  predictor-corrector,  but  both  of  these  are  difficult 
to  implement  with  white  noise  as  the  input.  Runge-Kutta  is  impossible  to 
implement  because  it  requires  the  evaluation  of  the  function  at 
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interscqucnce  values  of  time,  and  we  then  have  the  choice  of  using  a  new 
value  for  the  random  variable  or  of  using  the  value  at  one  end  of  the  time 
stop.  In  either  case  we  will  be  inputting  a  noise  spectrum  different  from 
that  of  the  discrete  white  noise  discussed  in  Section  3.2.  Predictor- 
corrector  methods  do  not  suffer  from  this  problem;  however,  when  they  are 
being  implemented,  the  values  of  the  noise  variable  at  the  several  time 
steps  involved  in  each  integration  step  must  be  consistently  maintained. 
This  involves  careful  programming,  so  it  is  preferable  to  use  a  simpler 
integration  routine  which  is  yet  well  behaved.  One  such  is  the  modified 
Euler  algorithm  given  by 


x(n  +  1 )  A 1  =  x[nA|  +  iA(x[(n  +  1)A]  +  x[nA]}  (31) 

This  is  also  known  as  second  order  Rungo-Kntta  algorithm.  It  gives  for  the 
discrete  transfer  function: 


A(1  +  z~M  _  _A  /j  +  costa  A 
2(1  -  z'1)  "  2 7  1  -  cosw  A 


(32) 


This  function,  which  is  also  plotted  in  figure  2,  is  seen  to  damp  out  high 
frequencies,  which  is  why  this  routine  is  stable.  Therefore,  provided  that 
the  integration  time  step  is  sufficiently  small  that  none  of  the  natural 
modes  of  the  system  is  severely  damped  by  the  transfer  function  of  this 
integration  routine,  it  is  a  suitable  integration  method.  It  also  has  the 
advantages  of  being  simple  to  implement  in  a  computer  program  (only 
requiring  two  iterations  per  step  interval)  and  of  being  simple  enough  to 
examine  analytically  to  see  what  effect  it  has  on  given  noise  inputs.  The 
well  known  tendency  of  integration  using  this  algorithm  to  'drift'  can  be 
seen  to  arise  because  the  curve  in  figure  2  starts  to  depart  from  unity  at 
the  lowest  frequencies. 

The  accuracy  of  the  modified  Euler  method  can  he  increased  by  iteration  for 
x(nA)  and  continuing  until  the  difference  between  successive  iterations  is 
as  small  as  desired.  When  this  is  done,  the  method  can  be  regarded  as  a 
second-order  prediction-corrector.  A  further  virtue  of  this  simple 
algorithm  is  that  it  is  suitable  for  a  set  of  stiff  equations(ref .6) ,  that 
is,  a  set  of  equations  with  a  wide  range  of  values  for  its  eigenvalues. 
Such  sets  of  equations  occur  in  flight  simulation  where  the  time  constant 
can  vary  from  100  s  (for  phugoid)  to  0.01  s  (for  actuators).  References  7 
and  8  discuss  the  benefits  of  the  modified  Euler  integration  routine  when 
the  set  of  equations  are  stiff,  it  is  frequently  superior  to  far  more 
sophisticated  algorithms. 

Although  it  is  not  the  purpose  of  this  memorandum  to  analyse  integration 
procedures,  we  will  end  this  Section  with  a  mention  of  an  alternative  to 
formal  integration  routines.  This  procedure,  discussed  in  reference  9, 
enables  large  time  steps  to  be  used  in  discrete  simulation  by  using  a 
difference  equation  which  accurately  maps  the  poles  of  the  original 
differential  equation.  This  technique  is  valuable  for  stiff  systems 
because  the  time  step  interval  can  be  made  comparable  with  the  time 
constant  of  the  fastest  mode  (instead  of  much  less,  as  is  required  for 
standard  integration  routines).  However,  when  the  coefficients  of  the 
difference  equation  are  derived  by  requiring ; correct  representation  of  the 
poles  of  the  original  differential  equation,  a  power  series  representation 
of  the  input  with  time  over  each  interval  must  be  assumed  (for  example, 
quadratic).  Since  white  noise  can  not  be  represented  by  a  power  series 
this  procedure  must  introduce  an  approximation.  It  is  generally  advisable 
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to  check  If  the  consequent  error  be  significant, 


Another  example  of  the  analogue-digital  technique  for  discrete  simulation 
of  continuous  systems  Is  given  in  reference  10.  In  the  method  given 
therein  the  coefficients  of  the  difference  equation  are  chosen  so  that  the 
first  and  second  moments  of  the  discrete  correlated  noise  are  the  same  as 
for  the  correlated  noise  generated  by  the  differential  equation  which  is 
being  modelled.  The  resulting  expressions,  which  are  rather  lengthy,  are 
essentially  in  terms  of  the  transition  matrix  of  the  differential  equation 
integrated  over  the  discrete  time  interval. 


We  will  now  study  the  generation  of  discrete  Markov  processes,  first  by 
difference  equations  derived  using  the  forward  difference  algorithm  and 
then  by  using  the  modified  Euler  algorithm. 

3. A  Generation  of  discrete  Markov  processes 


3.4.1  Forward  difference  algorithm 


By  using  the  forward  difference  approximation  to  the  derivative 
(equation  (29))  the  discrete  equivalent  of  the  first  order  Markov 
process  (equation  (12))  is 


x[(n+l)A|  =  (1-uA)  xfi'.Aj  +  nA.u(nA) 


(33) 


where,  as  has  been  shown  above,  u(nA)  is  a  sequence  of  uncorrelatcd 
random  numbers  with  variance  a 2  =  Y/A. 

The  power  spectral  density  of  x  is  given  by  the  product  of  the  modulus 
of  the  frequency  transfer  function  of  the  equation  with  the  power 
spectral  density  of  u(ref.ll),  so  we  have  for  the  discrete  first  order 
Markov  process  a  power  spectral  density: 


*XM 


_ a2  A2  * _ _ 

(z  -  1  +  vA)  (z-1  -  1  +  vA) 


(34) 


or 


*XW)  = 


_ a2  A2  * _ 

1  +  (1  -  vA)2  -  2(1  -  yA)  coscoA 


(35) 


Taking  the  limit  A-*0  gives  the  power  spectral  density  for  the  continuous 
first  order  Markov  process,  equation  (13),  as  expected. 

In  figure  3,  we  have  plotted  the  power  spectral  density  for  the 
continuous  Markov  process,  which  should  be  compared  with  that  of  the 
discrete  process  plotted  in  figure  4.  It  is  seen  that  the  power  for  the 
discrete  process  is  greater  than  for  the  continuous  process,  but  tends 
to  it  either  as  A^O  or  as  w*0. 


When  simulating  continuous  processes  by  discrete  sampling,  it  is 
important  to  know  what  errors  are  introduced  as  a  function  of  the 
sampling  interval.  We  will  now  derive  expressions  relating  the  maximum 
permissible  discrete  time  integration  interval  for  a  desired  accuracy 
-In- terms. of  tha. parameters  of  the  flrSt_order  Markov  process. _ ’ 
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There  are  two  constraints  on  the  sampling  interval,  A.  The  first  is 
that  the  power  spectral  density  decrease  to  a  negligible  value  at  the 
Nyquist  frequency  (=u/A) .  For  a  first  order  Markov  process,  the  ratio 
of  power  spectral  density  at  the  Nyquist  frequency  to  its  peak  value  is 

s  =  u2/ [  (tt/A)2+u2  ] 


We  can  interpret  this  equation  as  stating  that  a  fractional  error,  e, 
results  from  a  sampling  interval  A,  due  to  aliasing.  Rearranging  gives 


A  <  it/e/o 


(36) 


as  the  contraint  on  A,  where  we  have  used  the  condition  e<<1  in  deriving 
equation  (36).  The  other  constraint  is  that  the  sampling  interval  be 
smaxl  enough  that  the  discrete  power  spectral  density  approximates  the 
continuous  one.  Let  £  denote  the  ratio  of  the  difference  between  the 
discrete  and  continuous  pov  r  spectral  densities  to  the  maximum  of  the 
power  spectral  density.  This  ratio  maximises  at  w  =  u,  and  gives  the 
following  constraint  on  A,  for  th  error  ratio  to  be  less  than  e: 


A  <  2e/u 


(37) 


Note  that  the  two  constraints  both  involve  ».  This  is  hardly  surprising 
as  wo  have  only  one  parameter,  viz  u.  Since  £<<1,  equation  (37) 
imposes  tighter  constraint  or.  A  than  does  equation  t36). 

We  now  turn  to  second  order  Markov  processes.  Substituting 
equation  (33)  into  equations  (14)  and  (15),  eliminating  y,  and  taking 
the  2  transform  of  the  resulting  difference  equation  yields 


U1  +  2(i"to<>A-  l)z  +  (1  -  w§AJ)jx  =  (aA.r  *  (bA3  -  aA)]U  (38) 


where  b  -  c  +  2a4uQ  X,  U  arc  the  z  transforms  of  x  and  u  (the  z 

transform  of  u  being  purely  formal). 

We  can  now  write  the  power  spectral  density  of  r.  as 


4»  a 
X 


A3  [  { a'  »  ^bA  -  a)3j  ♦  J!a(bA  -  a)  costa  A] 

1  ♦  +  02  *  20i  fl  costa  A  +  202  cos2co  A 


(39) 


where 

0i  =  2(fw  oA-1) 

02  =  1  -  2$ to  0  A  +  to  q  A3 
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As  before,  wc  compare  the  continuous  second  order  Markov  process 
plotted  in  figure  6,  with  the  forward  difference  discrete  approximation’ 
f  I  ){iir<ts  7(a)  and  (b).  .It  is  seen  the  discrete  power  spectral  density 
diverges  from  the  continuous  more  rapidly  with  integration  step  interval 
than.it  did  for  the  first  order  Markov  case.  This  can  be  attributed  to 
a  A  resonance"  in  the  power  spectrum  (for  any  given  u) ,  regarding  $  now 
as  a  function  of  A,  and  is  illustrated  in  figure  7(c). 

As  before,  we  have  constraints  on  A  imposed  by  the  requirements  of  the 
continuous  power  spectral  density  being  small  at  the  Nyquist  frequency, 
and  of  the  discrete  power  spectral  density  being  close  to  the 
continuous.  These  are  more  difficult  to  apply  to  the  second  order 
Markov  process.  After  some  manipulation  it  can  be  shown  that,  for  the 
aliasing  and  discrete  approximation  errors  to  be  less  than  a  fraction  e, 
the  sampling  interval  must  satisfy 


A  <  jr  b</"c/aco  J 


(40) 


and 


A  < 


(b3  +  a3  to  p)f  e 
(b3  -  abfwo  +  a3  to  o  )to  0 


(41) 


3.4.2  Modified  F.uler  algorithm 

% 

We  will  now  derive  the  modified  Euler  approximations  to  the  first  and 
second  order  Markov  generating  equations  and  compare  them  with  the 
results  of  the  Section  3.4.1. 

Substituting  equation  (31)  into  equation  (12)  and  taking  the  z  transform 
gives 


[ (1+uV)  z  -  (1-uV)]  X  =  aV ( l+z)U 


(42) 


where  7  =  A/2. 

We  can  immediately  write  down  the  power  spectral  density  of  x  to  be 


*,w 


2a3  V3  fl  +  costoAl  ¥ 

1  +  vi  A3  -  (l  -  v2  A1)  coscoA 


(43) 


This  function  is  plotted  as  a  function  of  w  for  a  =  u  =  1  and  A  =  0,1,2 
in  figure  5.  Comparing  with  figure  4  we  see  there  is  a  substantial 
(about  four  times  in  this  case)  improvement  on  accuracy  for  a  given  A. 
Note  that  the  forward  difference  approximation  amplifies  the  power 
spectral  density  at  high  frequencies,  while  the  modified  Euler 
approximation  diminishes  it.  These  properties  are  to  be  expected  from 
their  respective  frequency  response  functions  (see  figure  2). 

Limitations  on  the  size  of  A  exist  from  the  same  constraints  as 
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discussed  in  Section  3.4.1  and  lead  to: 


A  <  We/ o 


A  <  3/e/u 


Comparing  equation  (45)  with  (37)  shows  the  modified  Euler  approximation 
allows  larger  time  intervals  for  a  given  accuracy. 

We  now  turn  to  our  final  case:  the  modified  Euler  approximation  to  the 
differential  equations  defining  the  second  order  Markov  process. 
Substituting  equation  (31)  into  equations  (14)  and  (15),  and  taking  the 
z  transform  gives 

1(1  +  2fw0V+  +  2tfu)  l  -  l)z  +  (1  -  2fw0V)]X 

=  V  [  (bV  +  a)  z2  +  2b V.  z  +  (bV  -  a)  ]  U  (46) 


whore  b  =  c  +  2aCw 

o 

which  gives  for  the  power  spectral  density 


a*  +  Sb^l  ♦  4b‘V*  cos  to  A  +  (b'y2  -  a2 )  cos  2oo  A]  ❖ 


03  *  04  +  0S  +  204(03  +  0s)  COS  t*i  A  +  2030s  cos  2toA 


who  re 


03  =  1  ♦  2f  to  oV  ♦  to  oV2 

04  s  2  (to  V2  -  1) 

0s  «  1  -  2f  to  oV  ♦  to  3v2 


Equation  (47)  is  plotted  as  a  function  of  u,  for  several  values  of  A,  in 
figure  8.  Comparing  with  figure  7  shows  that  the  modified  Euler  method 
is  able  to  use  much  larger  values  of  A  before  there  is  a  significant 
difference  between  the  continuous  and  discrete  spectra.  It  is 
interesting  that  the  modified  Euler  approximation  produces  not  only  a 
diminution  of  power  at  higher  frequencies,  but  also  a  lowering  of  the 
frequency  of  maximum  power.  In  place  of  equations  (40)  and  (41)  we  have 
the  following  inequalities  which  must  be  satisfied  for  the  discrete 
approximation  to  model  the  continuous  second  order  Markov  noise 
adequately: 


A  <  »b  y/H /&<■■'  2 
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6fb2  +  w  2  a2  )e 
.b1  (1  +  4?2)  ♦  d>  J  a1  . 


^  2fvST 


(49) 


We  have  summarised  the  results  of  Sections  2  and  3  in  Tables  1  through 
3.  We  have  included  the  expressions  for  .autocorrelation  function  in 
Table  1,  because  it  is  sometimes  this,  rather  than  the  power  spectral 
density  of  the  random  process,  which  is  measured  and  to  which  the 
theoretical  noise  process  is  fitted  to  determine  the  parameters  a,  u  or 
a,  c,  wq,  ?. 

3.5  Summary 

A  continuous  white  noise  stochastic  process,  with  power  spectral  density  ¥ 
(or,  equivalently  with  autocorrelation  function  ¥6^)  can  be  modelled  in  a 

discrete  simulation  by  a  sequence  of  uncorrelated  random  numbers  at  the 
simulation  times,  whose  variance  is  o*  =  ¥/A  where  A  is  the  discrete  time 
interval  of  the  simulation.  This  discrete  white  noise  also  has  a  uniform 
power  spectral  density  equal  to  ¥,  but  it  is  band  limited  to  the  angular 
frequency  range  -  n/A  <  w  <  it/A  (sec  figure  1). 

Although  white  noise  never  occurs  in  nature.,  its  benefit  lies  in  the  fact 
that  noise  of  any  given  power  spectral  density  can  be  generated 
analytically  by  the  solution  of  differential  equations  whose  forcing 
function  is  white  noise.  The  first  order  Markov  process  has  the  power 
spectral  density  (see  figure  3) 


<t>  = 


2  a2 

o2  +  cu 


2 


and  is  generated  by  the  first  order  differential  equation 


x  =  -ux  +  au 


where  u  is  the  white  noise  input,  and  is  always  chosen  so  that  ¥  =  1. 
The  second  order  Markov  process  (see  figure  6) 


a2 to  2  +  (c  +  2?  w  o  a)2 
x  =  cj4  +  2w  o  (2? 2  -  l)to2  +  lot  '  * 


is  generated  by  the  coupled  pair  of  equations 

x  -  y  +  au 

y  =  -u?  x  -  2f  «0y  ♦  cu 

Given  an  observed  power  spectrum  for  noise  we  are  then  able  to  model  the 
stochastic  process  which  has  this  power  spectrum  by  choosing  the 
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coefficients  of  the  equations  (viz:  a,  u;  or  a,  c,  w  ,  £)  so  the  resultant 
power  spectral  density  fits  the  observed  power  spectrum. 


In  discrete  simulation 
replaced  by  difference 
be  taken  in  choosing 
differential  operator, 
implementation  appears 


of  these  processes,  the  differential  equations  are 
equations.  As  is  discussed  in  the  text,  care  must 
a  suitable  algorithm  with  which  to  implement  the 
The  best  compromise  between  ■  "curacy  and  ease  of 
to  be  the  modified  Euler  method  r  which 


xl  (n+l)A]  =  xlnAj  +  JjAjxf  (n+l)Aj  +  xl  nA]  | 


where  the  derivatives,  x,  are  replaced  by  the  equations  to  be  integrated. 
The  numerical  procedure  which  implements  this  relation  is  a  two  step 
iteration,  because  an  estimate  of  x(n+l)  must  be  made.  When  writing  the 
code  for  the  discrete  simulation,  care  must  be  taken  to  be  consistent  in 
the  use  of  the  random  numbers  generated  to  simulate  the  white  noise  forcing 
function;  viz:  the  number  generated  for  u  at  the  second  step  in  the 

calculation  of  y(n+l)  should  also  be  used  in  the  first  step  in  the 

calculation  of  y(n+2). 

As  discussed  in  the  text,  integration  routines  such  as  fourth  order  Runge- 
Kutta,  which  require  calculations  at  times  intermediate  to  those  at  which 
the  difference  equation  is  evaluated,  arc  unsuitable  for  the  implementation 
of  noise.  This  is  because  of  the  problem  of  what  value  to  assign  to  the 

random  variable  at  the  intermediate  steps.  How  this  is  done  makes  a  large 

difference  to  the  power  spectral  density  of  the  process.  lienee  we 
recommend  simpler  integration  routines,  thereby  obviating  the  problem. 

A  further  matter  which  must  be  considered  is  that  the  power  spectral 
density  of  a  discrete  Markov  process  calculated  by  difference  equations 
differs  from  that  given  by  the  exact  solution  of  the  differential 
equations.  This  suggests  that  the  Markov  coefficients  should  be  calculated 
by  fitting  the  observed  spectrum  to  the  power  spectral  density  of  the 
discrete  process.  This  is  not  recommended,  because  the  fitment  is  no 
longer  valid  if  the  time  step.  A,  of  the  simulation  is  changed.  It  is 
preferable  to  fit  the  continuous  spectrum  to  the  observations  and  note  the 
error  introduced  by  discrete  simulation  for  various  time  stop  intervals. 
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TABLE  3.  DISCRETE  RANDOM  PROCESSES  -  GENERATION  BY  MODIFIED  EULER 
APPROXIMATION  TO  DIFFERENTIAL 
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NOTATION 

coefficient  of  white  noise  in  differential  equation  for  Markov 
process 

=  c  +  2  a  C  w 

o 

coefficient  of  white  noise  in  differential  equation  for  Markov 
process 

correlation  function 

expected  value  of  a  variation 

measured  value  of  a  state  variable  (=v+r*) 

continuous  noise  process 

probability  density  function 

pcv«r  (=  integrated  spectral  density  from  -»  to  «*) 

covariance  function 

Laplace  transform  variable  (=ju) 

time 

time  interval  for  integration  (in  the  limit  |T|'*«) 
white  noise  stochastic  process 
true  value  of  a  state  variable 

Markov  process  derived  from  white  noise 

(a)  first  order  =  exponentially  correlated  noise 

(b)  second  order  “  damped  cosine  correlated  noise 

intermediate  variable  in  deriving  second  order  Markov  process 
z  transform  variable  (=  exp(sA)) 

fb 

Dirac  delta  function,  f  6^(c-x)dx  s  1  if  a*c>l>;  -  0  otherwise 

J  a  U 

Kronecker  delta.  6^,^ f m)  s  1  if  a  a  0;  a  0  otherwise 
time  interval  of  discrete  process 
=  6/2 

fractional  error  in  discrete  approximation  to  continuous  power 
spectral  density 

damping  ratio  in  second  order  Markov  process 
mean  value 


correlation  time  in  first  order  Markov  process 
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o  standard  deviation  of  Gaussian  random  number  sequence  used  for 

discrete  simulation  of  white  noise 

t  time  lag,  argument  of  autocorrelation  function 

<l>  power  spectral  density  (two  sided  and  function  of  w) 

$  subsidiary  expressions  involving  £,  uq>  A 

w  Fourier  transform  variable.  Frequency  in  radian  s  1 

undamped  natural  frequency  in  second  order  Markov  process 

'I'  spectral  density  of  continuous  white  noise  (taken  =  1) 
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APPENDIX  I 

AN  EXAMPLE  OF  DISCRETE  NOISE  SIMULATION 

Suppose  wo  wish  to  model,  in  a  discrete  simulation  computer  program,  the  noise 
on  the  output  of  an  angle  measuring  device  whose  measured  noise  power  spectral 
density  is  given  in  figure  9.  This  approximates  to  a  second  order  Markov 
process,  and  so  we  will  assume  that  it  can  be  generated  from  white  noise,  u, 
by  the  following  pair  of  stochastic  differential  equations. 


Wo  x  -  0y+cu 

deg  s'2 

£1.1) 

y  +  au 

deg  s"1 

(1.2) 

As  shown  in  the  text,  the  power  spectrum  of  x  generated  by  these  equations  is 


$ 

x 


_ (a2  to2  +  b2) 

to  *  +  2to  o  -  l)to  2  +  to  q 


.  ❖ 


deg2 /(rad  s’*1) 


where  b  =  c  +  2a4WQ,  and  T  is  the  power  spectral  density  of  u.  Since  c  and  a 

appear  as  multipliers  of  u,  defining  ?  will  fix  c  and  a.  It  is  convenient  to 
take  t  as  unity,  and  we  will  therefore  do  this: 


viz: 


T  =  1  deg2  /(rad  s  l) 


The  parameters  a,  b,  u>o,  4  are  then  chosen  to  fit  4>  to  the  observed  power 

spectral  density.  This  is  partly  trial  and  error,  partly  intelligent 
guesswork;  for  example  the  frequency  and  relative  amplitude  of  the  maximum 
suggests  that  uq~2.5  and  4~ 0.5.  We  have  plotted  in  figure  9  the  calculated 

curve  for  the  following  values  of  the  parameters,  which  give  a  satisfactory 
fit : 


a 

= 

1.6  rad  s  1 

c 

= 

0,7  (rad  s  l) 

w 

0 

= 

2.6  (rad  s  *) 

4 

0.6 

b 

= 

c+2a4  w  - 
o 

To  generate  the  discrete  2nd  order  Markov  process  corresponding  to  x,  we 
replace  the  differentials  in  equations  ( 1 . 1 )  and  (1.2)  by  a  finite  difference 
algorithm.  As  explained  in  the  text,  a  good  compromise  between  convenience 
and  accuracy  is  given  by  the  modified  Euler  method.  Using  this  method,  and 
assuming  that  we  require  an  accuracy  of  better  than  20%  in  the  discrete 
modelling  of  the  continuous  process,  equations  (48)  and  (49),  give  the 
following  constraints  on  the  integration  step  size: 
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A  s.  ffl^e/au)  i  «  0.52  s 

A<2fvre/oo0  =  0.50  s 


We  have  plotted  in  figure  9  the  power  spectrum  of  the  discrete  process 
corresponding  to  a  0.5  s  time  step,  out  to  the  Nyquist  frequency 

(it/ A  =  6.3  sec-1)  . 

In  figure  10  we  have  plotted  100  s  of  the  discrete  random  process  generated  by 
the  modified  Euler  algorithm  applied  to  the  differential  equations,  with  the 
calculated  values  for  the  parameters,  generated  from  discrete  white  noise 
given  by  Gaussian  random  number  sequence  of  variance: 


c2  =  Y/A  =  2  deg2 


Tim  theoretical  variance  of  the  second  order  Markov  process  derived  by  these 
means  is  given  by  (see  Table  1): 


a2  _  (a2  co  o  +  b2)  y 

theory  4?w  o 

=  1.18  deg2 


This  value  should  be  compared  with  the  variance  of  the  calculated  random 
number  sequence  plotted  in  figure  10: 


a 


i 

calc 


1.0  deg2 


whore  the  sum  has  been  taken  over  500  steps.  It  is  wise  to  make  this 
comparison  as  a  chock  on  the  discrete  simulation  program.  These  values  could 
also  have  been  obtained  by  integrating  under  the  corresponding  curves  in 
figure  9.  The  20%  accuracy  between  the  continuous  and  discrete  power  is 
generally  adequate  for  most  applications.  It  is  not  often  that  the  estimate 
of  an  error  is  required  with  greater  accuracy. 

In  this  example  wo  started  with  a  given  power  spectral  density.  If  the 
measurements  of  the  noise  process  had  been  of  the  autocorrelation  function,  we 
could  determine  the  p> ^meters  a,  c,  C  by  fitting  the  theoretical  to 

observed  autocorrelation  function,  and  then  proceed  to  model  the  noise  in  the 
same  manner. 
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APPENDIX  II 

DISCRETE  MODELLING  OF  ATMOSPHERIC  TURBULENCE 

The  purpose  of  this  Appendix  is  to  demonstrate  how  the  methods  of  this  Report 
can  be  applied  to  the  simulation  of  atmospheric  turbulence.  In  simulating 
turbulence  it  is  necessary  to  use  an  analytical  expression  for  the  power 
spectral  density  of  the  lateral  turbulent  atmospheric  motions.  We  shall  apply 
the  Dryden  formula,  which  is  commonly  used;  it  has  the  advantages  of  being 
analytically  simple,  whilst  at  the  same  time  fitting  the  observed  spectra  very 
well. 

In  this  Appendix  we  do  not  intend  to  give,  nor  is  this  the  place  for,  an 
exegesis  on  atmospheric  turbulence,  its  variability  and  the  validity  of  the 
several  analytical  expressions  for  its  power  spectral  density.  These  topics 
are  covered,  for  example,  in  references  12  and  13. 

The  Dryden  expressions  for  the  power  spectral  density  of  the  longitudinal  and 
lateral  components  of  atmospheric  turbulence  are: 

"long  '  lVlhp  (■»-‘)1/(ra4.-)  (II. 1) 


1  +  3L2  ft2 

*ut  •  “  L  •  u  ♦  l!  a!P  (U-2J 


where 


Q  =  spatial  frequency  (rad  in"1) 

0*  =  mean  square  gust  velocity  (ms-1)2 

L  =  scale  of  turbulence  (m/rad) 


The  longitudinal  component  refers  to  those  random  atmospheric  motions  along 
the  line  to  which  the  spatial  frequency  is  referred;  the  lateral  component 
refers  to  motions  (vertical  or  horizontal)  perpendicular  to  this  line. 

If  the  above  equations  are  compared  with  the  Dryden  formulae  usually  quoted, 
it  will  be  noted  that  they  contain  an  extra  factor  of  ( 2v)/2 .  The  factor  2ir 

is  because  the  power  spectral  density  of  atmospheric  turbulence  is  commonly 

defined  so  that  its  integral  gives  the  variance,  o2.  However,  we  prefer  to 

retain  the  relation  between  power  spectral  density  and  variance  given  by 

equation  (7),  since  this  is  in  line  with  the  standard  definitions  of  the 
Laplace  and  Fourier  transforms.  The  factor  2  arises  because  we  are  using  a 
two-sided  power  spectral  dens' ty  of  to  be  consistent  with  the  main  text, 
rather  than  the  one-sided  density  which  is  commonly  used.  Finally,  note  that 
L  and  fl  are  in  radians  and  not  cycles. 

Numerous  observations  of  turbulence  have  shown  th«i  a  typically  lies  in  the 

range  0.5  to  3  ms"1,  with  the  lower  values  occurring  more  frequent ly(rcf .  14) . 
L  varies  from  150  m/rad  to  500  m/rad,  ex-  pt  near  the  ground  when  I.  is 
proportional  to  altitude.  Reference  12  recommends  using  300  m/rad  for  L  in 
the  free  atmosphere. 
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When  considering  the  effect  of  turbulence  on  an  aircraft,  it  is  standard 
practice  to  assume  that  the  aircraft  moves  through  a  spatial  field  given  by 
equation  (II.  1)  which  is  "frozen".  This  approximation,  called  the  Taylor 
hypothesis,  holds  as  long  as  the  aircraft  velocity  is  much  greater  than  all 
the  turbulent  motions.  The  aircraft  then  sees  these  spatial  structures  as  an 
equivalent  frequency,  u,  given  by 


u  =  VQ  rad  s"1 


where  V  is  the  velocity  of  the  aircraft.  Substituting  into  equations  (II.l) 
and  (II. 2)  gives  the  power  spectral  density  of  turbulence,  seen  by  an 

aircraft,  at  a  frequency  u  rad  s"1  as 


2p2L  1 

V  •  1  +  (Ia/V)5 


(ms"1)2 /(rad  s'1) 


(II. 3) 


o2 L  1  ♦  3(LoVV)2 
V  *  (  1  +  (Lu/V)2] 2 


(II. 4) 


We  now  wish  to  show  how  those  spectra  can  be  generated  by  using  the  methods 
developed  in  this  report.  Comparing  with  equations  (13)  and  (17)  it  is  seer, 
that  the  Dryden  spectra,  chosen  because  they  give  a  good  fit  to  observations 
while  at  the  same  time  being  similar  in  form  to  the  theoretically  exact  von 
Karmen  turbulent  spectrum,  can  be  modelled  exactly  by  the  first  and  second 
order  Markov  processes: 


for 

longitudinal  motions: 

x  = 

~ux  +  a’u  ms"1 

for 

lateral  motions: 

X  = 

y  +  au  ms"1 

y  a 

-oj  o x-2f<a  oy+cu  ms"2 

where  x  is  the  desired  output  and  u  is  white  noise  of  power  spectral  density 
? .  The  parameters  in  thoso  equations  are  related  to  those  of  the  Uryden 
spectra  by 


u  =  V/I. 

a”  =>  [  20JV/l)  /¥ 

C  =  1 

u  =  V/L 
o 

a1  =  (3o*V/Ll/T 
b1  »  (oaV*/L*J/T 
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where 


b  =  c  +  2a£w 

o 

and,  as  discussed  in  the  text,  we  let  the  white  noise  input  have  unit  power 
spectral  density  1  ¥  =  1  (ms  1)1/rad  s-1  . 

The  table  below  shows  the  value  of  these  parameters  for  L  =  300  m/rad  and 
three  values  of  RMS  gust  velocity,  as  observed  by  an  aircraft  travelling  at 

170  ms  l: 


0 

m  s"1 

0.5 

1.5 

3.0 

a' 

rad  s'1 

0.53 

1.60 

3.20 

a 

rad  s'1 

0.65 

1.96 

3.92 

b 

(rad  s"1)2 

0.21 

0.65 

1.29 

c 

(rad  s-1)2 

-0.53 

-1.59 

-3.19 

w  =  0.57  rad  s  1  for  all  o. 

o 


Note  that  this  approach  to  modelling  atmospheric  turbulence  produces  motions 
of  all  scales.  There  is  no  need  to  add  long  period  sinusoids  to  the  model  in 
order  to  simulate  the  large  scale  motions,  as  we  have  seen  done  elsewhere. 

For  completeness  we  will  note  the  autocorrelation  functions  corresponding  to 
the  Dryden  longitudinal  and  lateral  power  spectral  densities 
(equations  (II. 1),  (II. 2)}. 

Clong  =  0,o*PHC!/L)  (m/«)s 

Clat  =  °*(HCl/L)oxp(-UI/L)  (m/s)2 

These  equations  give  a  better  indication  of  the  meaning  of  L  than  do  the 
expressions  for  the  power  spectral  density. 
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APPENDIX  III 

DISCRETE  SIMULATION  OF  THE  WIENER  PROCESS 

The  Wiener  process,  also  known  as  random  walk  and  the  process  of  uncorrelated 
increments,  is  formally  defined  as  the  integral  of  white  noise  (references  1 
and  11): 

x  =  au  (III.  1 ) 


This  process  differs  from  those  examined  in  the  text  by  having  an 
autocorrelation  function  whose  expectation  does  not  converge  to  a  constant 
value  for  large  time.  On  the  contrary,  the  autocorrelation  of  the  Weiner 
process  increases  linearly  with  time: 


C ( t ,  t )  =  aaTt  t<t 


(HI. 2) 


where  t  is  the  elapsed  time  from  the  start  of  the  process,  and  we  will  again 
take  the  white  noise  process  to  be  of  unit  strength,  t  =  1. 

The  discrete  equivalent  of  the  Wiener  process,  using  the  modified  Euler 
integration  algorithm,  is 

x((n+l)A]  =  x(nA)  +  i  aA{u((n+l)A]  +  u(nA)}  (III. 3) 


where,  as  discussed  in  the  text,  u(nA)  is  a  sequence  of  random  numbers  whose 
variance  is  o1  *  t/A.  We  can  rewrite  this  equation  as 


xl  (n*l)Al  =  aA)  u(o)  ♦  )  u(nA)  ♦  ul  (n+l)A] 


The  autocorrelation  function  can  therefore  be  written  as,  letting  m<n: 


C^Cra.n)  =  K(x(mA)  x(nA)j 

=  EE(aA)1  E(u(mA)  u(nA)) 
3  (aA)1  IZ  os  -  n) 

=  (m  -  maoA)J 


where  the  *  arises  from  the  end  point  weighting  of  u. 

Substituting  the  relation  o*A  =  t  and  setting  t=«A  gives 

C^(ra,n)  =  (m-$)/o.fiJYt  m<n. 

We  see  that  discrete  simulation  gives  the  same  autocorrclot ion  as  the 


continuous  process,  except  for  the  factor  (m  -  })/m,  ie  the  accuracy  of  the 
discrete  simulation  increases  with  the  number  of  steps.  The  step  size  can 
therefore  be  chosen  to  give  the  desired  accuracy  at  any  specified  time. 

Because  of  the  time  dependent  property  of  the  autocorrelation  function,  the 
power  spectral  density  is  not  defined. 
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Figure  2 


Hie  curves  give  the  ratio  of  the  frequency  response  of  the 
given  discrete  integration  algorithm  to  the  frequency  res¬ 
ponse  of  exact  integration  (=  1/cJ) 


Figure  2.  Comparison  of  frequency  response  of  two  discrete  integration 
algorithms 


"El 


d  order  Msrkov  process 


(c)  Variation  with  step  size 


Figure  7 


Power  spectral  density  of  second  order  Markov  process:  Discrete 
with  forward  difference  algorithm 


§  (°)2 /  (RAO 
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Figure  9 


_ - _ measured  power  spectral  density  (two  sided) 

—  — *  fit  of  2nd  order  continuous  Markov  process  to  measured 
p.s.d.  (see  Appendix  I) 

. -  discrete  approximation  to  —  with  A  =  0.5 


Li  (s) 


Figure  9.  Measured  power  spectral  density  fitted  by  second  order 
Markov  process 
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